Multiplex immunofluorescence and single‐cell transcriptomic profiling reveal the spatial cell interaction networks in the non‐small cell lung cancer microenvironment

Abstract Background Conventional immunohistochemistry technologies were limited by the inability to simultaneously detect multiple markers and the lack of identifying spatial relationships among cells, hindering understanding of the biological processes in cancer immunology. Methods Tissue slices of primary tumours from 553 IA∼IIIB non‐small cell lung cancer (NSCLC) cases were stained by multiplex immunofluorescence (mIF) assay for 10 markers, including CD4, CD38, CD20, FOXP3, CD66b, CD8, CD68, PD‐L1, CD133 and CD163, evaluating the amounts of 26 phenotypes of cells in tumour nest and tumour stroma. StarDist depth learning model was utilised to determine the spatial location of cells based on mIF graphs. Single‐cell RNA sequencing (scRNA‐seq) on four primary NSCLC cases was conducted to investigate the putative cell interaction networks. Results Spatial proximity among CD20+ B cells, CD4+ T cells and CD38+ T cells (r 2 = 0.41) was observed, whereas the distribution of regulatory T cells was associated with decreased infiltration levels of CD20+ B cells and CD38+ T cells (r 2 = −0.45). Univariate Cox analyses identified closer proximity between CD8+ T cells predicted longer disease‐free survival (DFS). In contrast, closer proximity between CD133+ cancer stem cells (CSCs), longer distances between CD4+ T cells and CD20+ B cells, CD4+ T cells and neutrophils, and CD20+ B cells and neutrophils were correlated with dismal DFS. Data from scRNA‐seq further showed that spatially adjacent N1‐like neutrophils could boost the proliferation and activation of T and B lymphocytes, whereas spatially neighbouring M2‐like macrophages showed negative effects. An immune‐related risk score (IRRS) system aggregating robust quantitative and spatial prognosticators showed that high‐IRRS patients had significantly worse DFS than low‐IRRS ones (HR 2.72, 95% CI 1.87–3.94, p < .001). Conclusions We developed a framework to analyse the cell interaction networks in tumour microenvironment, revealing the spatial architecture and intricate interplays between immune and tumour cells.

• Deep learning algorithm on multiplex immunofluorescence images identified cell spatial patterns with prognostic effects in the lung cancer microenvironment.
• Single-cell RNA-sequencing revealed the cell interaction networks suggested by spatial paradigm analyses.
• Proximity among CD4+ T cells, CD20+ B cells and N1-like neutrophils, and spatial compartmentalisation between cancer stem cells and CD8+ T cells, were associated with significantly longer disease-free survival.

INTRODUCTION
The tumour microenvironment (TME) contains both cancer cells and multifarious lineages of stromal cells, including various assemblages of immune, vascular and mesenchymal cells, contributing to the de novo development of cancers. 1,2 Emerging evidence showed that, apart from the amounts and types of cells within the TME, the spatial configuration is crucial for modulating response to therapies and survival in various types of tumours, such as lung cancer (LC), 3 breast cancer (BC) 4 and gastric cancer. 5,6 Approaches to comprehending the intricate biology and the predictive and prognostic connotations of different cell populations in TME are mushrooming. Research focuses are transitioning from single-color immunohistochemistry (IHC) imaging assays to multiplex IHC and multiplex immunofluorescence (mIF) techniques, enabling simultaneous detection of diverse protein markers in one tissue slice for the identification of multiple cell phenotypes. 7,8 Given that conventional IHC assays are confronted with several limitations like the inability to simultaneously detect multiple markers per tissue slice and the lack of identification of spatial relationships among cells, methods incorporating spatial organisation analyses are further developed. For instance, the metrics which quantify this spatial architecture, like nearest neighbour distance 9 and mixing scores, 4 have been proposed recently. Understanding the location of cells facilitates deeper insights into how immune cells interplay with tumour cells because both indirect and direct cell signal mechanisms require proximity between cells. 10 LC is the most common cause of cancer death globally. 11 Despite significant progress in new medicines and systemic therapies, like immune checkpoint inhibitors targeting TME, patients' overall survival rate was still dismal. 12 Consequently, there is an urgent need to improve the understanding of the TME of LC to better stratify patients for treatments and recognise novel targets for therapies to improve outcomes. The advances in spatial paradigm analyses, like deep convolutional neural networks, allow for a better description of cellular proximity in TME. 13 Wang et al. 14 developed a prognostic model of lung adenocarcinoma (LUAD) based on spatial features extracted by a deep learning algorithm on hematoxylin and eosin (HE)stained pathological graphs, improving the understanding of spatial locations of different cells and patient outcomes. However, they did not categorise the identified lymphocytes into different sub-populations, hampering the interpretation of the interactions among different lineages of lymphocytes. A study by Enfield and colleagues 3 reported that CD3+ CD8− T cells and CD8+ T cells surrounded by tumour cells were correlated with no recrudescence of LC, underscoring the cell sociology patterns play a crucial part in anti-tumour immune response. Nonetheless, the sample size of their study was rather limited (n = 20), which may lack generalisation of the results.
Herein, we investigated the spatial immune-tumour cell communication patterns and prognostic implications through stepwise procedures. First, the mIF test detecting 10 markers from 553 non-small cell lung cancer (NSCLC) cases was conducted, evaluating the amounts of 26 phenotypes of cells in tumour stroma (TS) and tumour nest (TN).
Second, the StarDist depth learning model was trained and developed using more than 2700 whole-slide mIF images, generating the features to depict the TME from the identified phenotypes and spatial location of these cells. 15 Eventually, to further understand the potential biological processes underlying the prognostic effects of different spatial patterns, we profile the transcriptomic features of the primary tumours and explore the putative cell interaction networks through receptor-ligand and downstream function analyses based on single-cell RNA sequencing (scRNA-seq). We found complex interactions among immune and tumour cells and revealed novel spatial parameters that were associated with the recurrence of LC. We further developed a prognostic tool, the immunerelated risk score (IRRS), to predict disease-free survival (DFS) across different stages, resulting in stable and great performance.

Study population and tissue collection
Formalin-fixed, paraffin-embedded (FFPE) tumour samples of NSCLC were retrospectively collected from 553 patients who underwent sub-lobectomy/lobectomy with lymphadenectomy in the First Affiliated Hospital of Guangzhou Medical University (GZMU) between 2009 and 2011.
Patients inclusion criteria were: (a) staged IA to IIIB primary NSCLC based on the National Comprehensive Cancer Network (NCCN) criteria 16,17 ; (b) resected tissues were sufficient for mIF test; (c) pathological confirmation of all lymph nodes and resected tissues were provided; (d) without macroscopically or microscopically positive surgical margins. Patients exclusion criteria were: (a) patients with non-invasive LC (e.g., minimally/ in situ invasive LUAD); (b) patients with preoperative neoadjuvant therapy; and (c) preoperative diagnostic biopsy.
The period from radical excision to local recurrence was defined as DFS. The ethics committee of First Affiliated Hospital of GZMU approved this study, and the study was carried out based on the Declaration of Helsinki. 18 All patients provided informed consent for sample collection and mIF analysis, and all participants consented to the publication of research results.

Multiplex immunofluorescence test
The mIF staining was completed at Genecast Biotechnology Co., Ltd. (Beijing, China). Two panels, including a total of 10 markers, were detected, which were CD38, CD20, CD4, FOXP3, and CD66b in panel 1, and PD-L1, CD163, CD8, CD68, and CD133 in panel 2. Multiple images were from serial sections from the same block per patient, and they were stained by DPAI and five markers in panel 1 or panel 2. Detection of each panel was based on a 4-μm thick slide cutting from FFPE NSCLC tissues. After deparaffinisation and rehydration, the slides were subjected to epitope retrieval through boiling for 20 min in Tris-EDTA buffer (pH = 9; Klinipath #643901, the Netherlands) at 97 • C. Subsequently, endogenous peroxidase was blocked by incubation for 10 min in Antibody Block/Diluent (PerkinElmer #72424205, USA), and later the protein was blocked in 0.05% Tween solution at 26 • C for 30 min. Then the five antigens in each panel were labelled by cyclic staining, including incubation of primary and secondary antibodies, tyramine signal amplification (TSA) visualisation and removing the TSA-antibody complex in Tris-EDTA buffer by microwave treatment for 20 min at 97 • C. In each round, antibody labelling was followed after epitope retrieval and protein blocking as mentioned above. After cycle staining, each slide was counterstained with DAPI for 5 min and mounted in Pro-Long Diamond Antifade Mountant (Thermo Fisher). As a positive control, fresh whole-slide of normal human tonsils were included in every staining batch to evaluate the experiments' reproducibility.
Furthermore, primary antibodies for CD133 and CD4 were incubated overnight at 4 • C; CD20, CD38, PD-L1, CD66b, CD8, CD68, CD163 and FOXP3 were incubated for 1 h at indoor temperature. Information of primary antibodies used was provided in Table S1. For secondary antibodies, we utilised anti-rabbit/mouse horseradish peroxidase antibodies (Zsbio # PV-6002 or PV-8000) and incubated them for 10 min at 37 • C. TSA visualisation was carried out through the Opal seven-colour mIF Kit as we previously described. 19 All slices were scanned utilising the PerkinElmer Vectra (Vectra 3.0.5). First, we used a 4× microscope objective lens to preview the whole picture of the slide, and then the 20× microscope objective lens was utilised to capture the details of the visual field. Finally, the images of each field were spliced to obtain the full view of the slide. The images were of high resolution (4028 × 3012 px), contributing to precise detection. Through the inform Advanced Image Analysis software (version 2.3.0), multi-spectral images were unmixed with spectral libraries constructed from images with single stained tissue for each antigen.
First, 10-15 mIF images with high resolution in the image dataset were selected by random sampling (Supplementary Data 1). Second, our experienced pathologist (Dr. Bai Xuejuan) delineated the areas of TS and TN on the images to train the algorithm in the inForm software. Third, the inForm software was capable of detecting and segmenting the tissue slides into TS and TN according to their morphologies automatically. Segmentation of cells was also carried out in the inForm software ( Figure S1). Determination of the suitable positive threshold X of each marker was implemented by two experienced pathologists (Dr. Bai Xuejuan and Dr. Wang Xin) independently, and disagreements were settled by consensus. Respectively, X, 2X and 3X was defined as the threshold of low ('+'), median ('++') and high fluorescence strength ('+++'). As illustrated by our previously article, 20 the histochemistry score (H-score) was conducted as: H − score = (cells with low f luorescence strength) % × 1 + (cells with median f luorescence strength) % × 2 + (cells with high f luorescence strength) % × 3 Markers for annotation of different cell were available in Table S2. 21

Identification, segmentation and localisation of cells by deep learning model
The dataset included more than 2700 whole-slide pathological images of 553 LC patients, and multiple images for one patient were from the serial sections from the same block. The deep learning pipeline to identify, segment and locate cells is illustrated in Figure 1.
We attempted to extract the fluorescent and get the cell coordinates after removing the noise because some nontarget cells may also be stained. Our fluorescent staining images contained R, G and B (i.e., red, green and blue, respectively) three channels. Given that the recognition performance was better in the R and G channels than in the B channel, then we separated the G and R channels of the fluorescent staining images, identifying six (CD4, CD20, CD38, CD66b, CD133 and CD163) and four (CD8, CD68, PD-L1 and FOXP3) types of markers, respectively. Finally, we extracted the fluorescently stained areas for each image. For DAPI-stained images, we used the StarDist depth training model to identify, segment and locate nuclei. 15 The model predicted the polygon corresponding to each pixel, then parameterised the radial distance to obtain the target probability. Then U-Net network was used to predict convex polygons, and finally, a final set of polygons was selected by non-maxima suppression, each polygon representing a separate object instance. We inputted the DAPI-stained image into the StarDist model to obtain the coordinates of the nucleus, that is, the cell coordinates, and then take the intersection of these cell coordinates and the coordinates of the fluorescent staining region extracted in the previous step to obtain the specific coordinates of each fluorescently stained cell after denoising.
We use the Delaunay triangulation algorithm to connect cells. The feature of this algorithm is to ensure that the nearest points form a triangle, the sum of the side lengths of the triangles is as small as possible, and each Delaunay triangle is circumscribed. The circle does not contain any other points on the surface. It also ensures that the minimum interior angle of the triangle in the Delaunay triangulation that the point set can form is as large as possible. And the triangle is as close to an equilateral triangle as possible. The connection diagram obtained by this connection method is reasonable and is more in line with the biological characteristics of cells in the complex TME.

Evaluation of segmentation performance
The ratio between the identified nuclei and the total ground truth nuclei was defined as detection coverage. Every ground truth nucleus matched the segmented nucleus, generating the maximum intersection over union (IoU). It is acknowledged that, at the object level, the higher the F1-score is, the higher the model's segmentation quality. At the same time, the IoU is a crucial factor contributing to the F1-score. Schmidt et al. 15 have developed the Stardist model based on the pathological slides, and it achieved optimal cell segmentation performance with the threshold value of IoU set as 0.6. Given the large sample size (ours vs. Schmidt U's: 2700 images vs. 1000 images) and high resolution (ours vs. Schmidt U's: 4028 × 3012 px vs. 256 × 256 px) of our image dataset, it was time consuming and inefficient to label all the ground truth nuclei. Consequently, utilising the same Stardist model in the present study, the threshold of IoU was also set as 0.6. The nucleus was marked 'matched' if IoU for a ground truth nucleus > 0.6; otherwise, it was marked 'unmatched'.
The pixel accuracy (PA) and F1-score were used to assess the segmentation accuracy. Two experienced pathologists (Dr. Bai Xuejuan and Dr. Wang Xin) manually labelled the cells on 15 pathological images with random sampling and localised the nuclei by star-convex polygons in the Stardist model. Combining the labels of cells provided by pathologists and the results predicted by the algorithm, we counted the values of true positive (TP), false positive (FP) and false negative (FN) ( Figure S2). The accuracy and recall of cell nucleus segmentation were also calculated to generate the F1-score based on the formula 22

Features extraction to describe the spatial relationships among cells
We extracted the number of connections and spatial distances between different types of cells. According to the length of each type of cell list, the subscript intervals of the coordinates of different types of cells in the general table could be obtained. The merged coordinates were inputted into the triangulation algorithm for connection, and the list of each edge in the output triangulation is called an edge table. Each edge in the edge table was a small list consisting of the numbers of the two endpoints that the edge corresponds to.
Here, we purposely matched the number of the input coordinates of the triangulation algorithm with the subscript of each element in the general table. Each edge in the table corresponded to the cell type and specific coordinates of the endpoint. On this basis, we established a loop to traverse all the edges on the triangulation, quickly determining the cell type and coordinates of the two endpoints of each edge and performing connection statistics and Euclidean distance calculations. After the loop was over, the number of connections between the two cells was output, and the spatial distance was obtained by dividing the total distance of each connection between the two cells by the number of connections.
Because 10 markers, including CD66b, CD4, CD38, FOXP3 and CD20 in panel 1, and CD163, PD-L1, CD8, CD133 and CD68 in panel 2, were enrolled in the present work, the edges of the image were classified into 30 categories based on the formula of combination: The number of connections (i.e., edges) for different classifications (yielding 30 features) and the average lengths of the connections (another 30 features) for each edge category were calculated for each graph patch. Eventually, 60 image features were extracted in total. The amounts of features were averaged for each patient if more than one pathology slide of one patient were available.

Quantitative and spatial features of cells in tumor microenvironment
Immune infiltration profiles, including the amounts and spatial distribution of cells, were used to decipher the TME characterisation. The percentage (%/sight) and density (per mm 2 ) of cells in both TN and TS were measured by the mIF test, while the spatial variables, including the number of connections and the average length of connections between two types of cells ( Figure S3), were summarised by mIF image-derived deep learning model. Noteworthily, given that both indirect and direct cell signal mechanisms require proximity between cells, more connections and shorter length of connections between two types of cells indicate closer interplays between them. In contrast, fewer connections and longer lengths of connections suggest weakened interactions. We further calculated the relative amounts of spatial variables in each mIF slice for downstream analyses: A total of 66 cell amounts variables and 60 spatial variables were yielded ultimately. Through the Spearman rank correlation test, correlations between two variables were evaluated to explore the distribution patterns and possible communications.

2.7
Generation of the immune-related risk score model The entire cohort (n = 553) was split into a training cohort (n = 387, 70%) and a testing cohort (n = 166, 30%) using stratified sampling referring to the 10 events per variable (10 EPV) principle. 23 Using the univariate and multi-variate Cox regression analyses, markers significantly associated with DFS were identified. Second, a Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression model was fitted to select the robust prognosticators among candidate factors with prognos-tic ability in univariate Cox results. Utilising the 10-fold cross-validations, we obtained the optimal lambda with the minimal partial likelihood deviance value. 24 Finally, the proportion of the selected prognosticators and their regression coefficients of the univariate Cox regression were multiplied to construct the IRRS model. The IRRS for each patient was calculated as follows 25 : wherein HR z is the HR for prognosticators, and proportion z represents the proportion for the z-th variables. To evaluate the prognostic performance of the IRRS system, multi-variate Cox regression analysis was applied. The assignment of high-IRRS and low-IRRS groups within the training cohort was executed following the optimal cut-off value. In the testing and entire cohort, identical formulas and cut-off value were utilised to validate the robustness of IRRS model.

Implementation of the single-cell RNA sequencing
Four primary invasive adenocarcinoma samples from four patients after surgical resection (Table S3) were rapidly transported in ice-cold H1640 (Gibco, Life Technologies) culture media. The tumour tissues were rinsed with phosphate-buffered saline (PBS), cut into 1 mm 3 cubic pieces, and digested with 0.25% trypsin. After being terminated by culture media supplemented with 10% foetal bovine serum, dispase (0.6 U/ml) and 10 ml of digestion medium with collagenase IV (100 U/ml) were added to the specimens. The digested samples were filtered by a 70μm nylon mesh, and filtered cells were collected through centrifuging for 5 min at 120×g, and 4 • C. Cell pellets were re-suspended with ice-cold red blood cell lysis buffer and subjected to a 40-μm nylon mesh. Following centrifugation, the cells were then collected with PBS, and the number of live cells were then calculated by an automatic cell counter. Cells were maintained on the ice during the process, and the entire process was completed in less than 1 h. Single-cell suspensions were switched to barcoded scRNA-seq libraries by the Gel Bead Kit, Single Cell 3′ Library, and Chromium Single Cell A Chip Kit (10X Genomics). To generate single-cell gel beads in the emulsion (GEMs), the cell suspension was recorded into the Chromium single-cell controller. The released RNA from lysed cells was barcoded by reverse transcription reaction in individual GEMs, and complementary DNA was then generated. The scRNA-seq libraries were then established utilising the Single Cell 3′ Library Gel Bead Kit V3, and sequencing was done by an Illumina NovaSeq 6000 with a paired-end 150-base pair (PE150) reading strategy (CapitalBio Technology, Beijing).

2.9
Processing single-cell RNA sequencing data Using human reference version GRCh38, raw expression matrices of each sample were created by the Cell Ranger (version 6.0.1) pipeline, and the Seurat package (version 4.0.3) was employed to analyse output-filtered gene expression matrices. 26 Once the following criteria were met, the cells were classified as low-quality ones: (i) <200 or >5000 genes, (ii) <500 unique molecular identifiers (UMI) or (iii) >10% UMIs derived from the genome of mitochondria. After removing the low-quality cells, the NormalizeData function was utilised for the normalisation of the gene expression matrices, and the FindVariableFeatures function was applied for calculating 2000 features with great intercellular variability. To reduce the dimensionality of the datasets after linear transformation with default parameters by the ScaleData function to generate scaled data, the RunPCA function was employed. Finally, all cells were clustered by the FindClusters and the FindNeighbors functions, and non-linear dimensional reduction was conducted through the RunUMAP function. In brief, after the identification of 2000 features with great intercellular variability, the FindIntegrationAnchors function was used to generate the anchors between individual datasets. The IntegrateData function was then applied with the anchors to produce an expression matrix of cells with batch-corrected from different datasets used for further analyses.

2.10
Annotating cell types and identifying cluster markers Through a uniform manifold approximation and projection (UMAP) approach after nonlinear dimensionality reduction, cells from primary tumours were projected into the two-dimensional field. Cells sharing similar characteristics could be clustered together subsequently. Widely accepted and classical markers of the corresponding cell types were utilised to annotate the clusters manually. 26

Profiling interaction networks among different cell types
The receptor-ligand analysis based on the CellChat algorithm, which could evaluate the intercellular interaction networks from scRNA-seq data, was implemented to assess the putative communications among different cell types. 27 The Secreted Signaling Pathways and the Cell-Cell Contact databases were selected for the main analyses. Two vital functions, including computeCommunProb and netVisual_bubble, were carried out to compare the probabilities of interactions mediated by ligand-receptor pairs from a certain cell group to another. To validate the observed communication networks, we obtained external processed scRNA-seq data of LUAD (n = 13) in the GSE123904 cohort from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and LUAS (n = 2) in the E-MTAB-6149 cohort from the Array-Express database (https://www.ebi.ac.uk/arrayexpress/) for the same analyses. 28,29

Gene set variation and enrichment analysis
To investigate the differences of function and allocate differential activities of pathways between different cell lineages, gene set variation analysis (GSVA) (R version 1.40.1) in hallmark gene sets and gene set enrichment analysis (GSEA) (version 1.44.0) in gene ontology sets (GO) were carried out. [30][31][32] Pathways with a p_adj value less than 0.05 were considered significantly enriched.

Statistical analysis
Spearman rank correlation test was exploited to explore the associations between spatial and amounts variables.

Demographic characteristics of included patients
A total of 553 patients with stage IA-IIIB NSCLC that met inclusion criteria were eventually enrolled, among which 240 cases (43.4%) relapsed during the follow-up time ( Table 1). The median DFS was 1678 days, and the 1-, 3-and 5-year recurrence rates were 13.9, 36.9 and 60.4%, respectively. Male (58.8%), stage IB (30.2%) and LUAD (69.1%) patients constituted the majority of the cases. Patients with more advanced clinical (cTNM) stage (p < .001) and confirmed visceral pleural invasion (VPI) (p = .015) had a greater propensity for recurrence (Table 1).

Segmentation performance of the StarDist model
In the pixel-level accuracy, the value of PA was 99.996%, and the F1-score was 81.39% in the object-level accuracy. In conclusion, using the StarDist model to identify, segment and locate the cells was effective and stable.

Deciphering immune landscape and interactions among cells
It is acknowledged that cell-cell interaction is an integral part of the biological process of antitumour/protumour immunity, exemplified by the intricate network of communications between target cells and diverse immune cells. We first explored the correlations between different cell populations by their amounts. We found that CD38+ T cells and CD133+ cancer stem cells (CSCs) were the most abundant cell populations in TS and TN, respectively (Figures 2A and B).  Figure S4), similar to our previous findings. 19 With the development of spatial transcriptome study, it is anticipated that the spatial analysis can produce more useful and advanced information which cannot be answered by rough measures of cell amounts. Consequently, we then investigated the cell distribution patterns to better decipher the TME (Table S4). Closer proximity and more connections were observed in CD20+ B cells, CD38+ T cells and CD133+ CSCs. In comparison, longer distances and fewer connections were found in CD68+ macrophages and CD8+ T cells, and neutrophils and CD4+ T cells (Figures 2C-E).
It is well known that regulatory T cells (Tregs), marked by FOXP3, are the major contributors to forming immunosuppressive TME by releasing cytokines like transforming growth factor (TGF)-β and interleukin (IL)-4. 33 Accordingly, we found modest and inverse associations between the distribution of FOXP3+ Tregs and CD20+ B cells. That is, the closer proximity between Tregs predicted fewer CD20-CD20 connections (r 2 = 0.32) ( Figure 2F). And fewer connections between CD20+ B cells indicated weaker interactions, which may lead to their dysfunction. The distribution of FOXP3+ Tregs was additionally correlated with decreased infiltration level of CD38+ T cells (a subset of CD4+ T cells) (r 2 = −0.36) ( Figure 2G). 34 Interestingly, proximity between CD66b+ neutrophils and FOXP3+ Tregs were associated decreased infiltrating levels of CD20+ B cells and CD38+ T cells (r 2 = −0.45) ( Figure 2H). Given CD20+ B cells and CD38+ T cells are involved in activating cytotoxic T cells for killing tumour cells, such distributed pattern seemed to show that FOXP3+ Tregs and CD66b+ neutrophils populations may synergistically attenuate this process (average length of CD38-FOXP3 connections and number of CD20-CD20 connections, r 2 = 0.39; number of CD38-CD66b connections and number of CD20-CD20 connections, r 2 = −0.31). Moreover, we found proximity among CD4+ T cells, CD20+ B cells and CD38+ T cells, and such spatial structure has been reported to aid the formation of tertiary lymphoid structure (TLS) 35 (Figures 2I-K).
Tumour-associated macrophages (TAMs), mainly M2 macrophages (marked by CD163), are pivotal players in shaping TME. 36 It has been reported that CSCs can produce TGF-β in the surrounding milieu, favouring the functional polarisation of protumour TAMs. 37 And in turn, TAMs facilitate protecting CSCs from the assaults of immune cells through expressing higher levels of PD-L1 and releasing cytokines like IL-10 and TGF-β. 38 Accordingly, we found close and positive associations between the distribution of macrophages, CD133+ CSCs and PD-L1+ cells (average length of CD133-CD163 connections and average length of CD68-PDL1, r 2 = 0.32) ( Figure 2N). Additionally, we found proximity among macrophages, CD133+ CSC and CD8+ T cells (average length of CD8-CD133 and average length of CD163-CD163, r 2 = 0.39; average length of CD133-CD163 connections and average length of CD8-CD68 connections, r 2 = 0.41; average length of CD8-CD163 and average length of CD68-PDL1, r 2 = 0.32) (Figures 2L-O), and such spatial pattern has been shown to prime T cells for exhaustion, possibly by the longterm synaptic interactions between them. 39 Ultimately, consistency analysis showed that higher infiltration levels of CD8+ T cells or CD133+ CSCs were associated with closer proximity between them (Figures S5A-E). And such spatially distributed patterns were supported by previous studies that CSCs tend to distribute closely to form the 'niches', which facilitate protection against the assault of the immune system and reserve their potential to metastasize. 40

Differences in cell infiltration and spatial location between lung adenocarcinoma and squamous carcinoma
LUAD showed significantly higher infiltrating levels of CD4+ T cells in TS, while lung squamous cell carcinoma (LUSC) was characterised by higher infiltrating levels of PD-L1+ cells, CD133+ CSCs and macrophages both in TS and TN, suggesting immune suppressive milieu ( Figure  S6). Closer proximity and stronger interaction between CD8+ T cells were observed in LUAD, whereas LUSC featured spatial proximity among PD-L1+ cells, CD133+ CSCs and macrophages. In conclusion, T cells, CSCs and macrophages were crucial contributors to the heterogeneity in TME between LUAD and LUSC, aligning with the recent findings by Wang et al. 41

Differences of cell spatial distribution patterns between patients with or without recurrence
The proximity between CD133+ CSCs was observed in patients with tumour relapse than without ( Figures 3A  and B). Concurrently, longer distances between CD66b+ neutrophils ( Figure 3C), CD66b+ neutrophils and CD4+ T cells ( Figure 3D), CD66b+ neutrophils and CD38+ T cells ( Figure 3E) and CD66b+ neutrophils and CD20+ B cells ( Figure 3F) were observed in patients with recurrence than without, suggesting the putative vital effects of neutrophils on spatially adjacent T and B lymphocytes. Apart from these cell types, different spatial landscapes in patients with or without relapse were observed ( Figures 4A and B).

Associations between spatial features and clinical characteristics
The distance between CD20+ B cells and CD66b+ neutrophils gradually increased with the advancing T stage (Figures S7A-F). Although null prognostic significance was found by multi-variate Cox regression analysis, the distance between CD38+ T cells and CD4+T cells increased along with the more advanced cTNM and N stage ( Figure  S7D).

Construction and validation of the immune-related risk score model
A total of 35 variables, including 30 cell amounts variables ( Figure S8C) and five spatial variables, which were strongly predictive of DFS in the univariate Cox analyses, were incorporated into the LASSO regression model and ten-fold cross-validation to generate the best model. Eventually, the IRRS model incorporating five amounts variables and five spatial variables achieved the optimal efficiency to speculate the DFS (Figures 5A and B).
The identified robust prognosticators could be crudely classified into several categories: antigen processing and presentation (intratumoural CD68-positive macrophages, average length of CD4-CD20 connections, average length of CD20-CD66b connections),

F I G U R E 3
Comparison of spatial cell location in non-small cell lung cancer (NSCLC) cases with or without recurrence. Discrepancies of the spatial distribution of cells between recurrent and non-recurrent patients as evaluated by the Kruskal-Wallis H test (A). Representative multiplex immunofluorescence images exhibited the spatial configuration of cells within two fields from one tissue slide. More number of CD133-CD133 connections (B), longer distances of CD66b-CD66b connections (C), longer distances of CD4-CD66b connections (D), longer distances of CD38-CD66b connections (E) and longer distances of CD20-CD66b connections (F) were significantly associated with NSCLC recurrence. *p < .05; **p < .01; ***p < .001. H test (A and B). The prognostic effects of spatial features, including the number of connections (C) and the average length of connections (D), in disease-free survival (DFS) as evaluated by univariate Cox analysis. *p < .05; **p < .01; ***p < .001; ns, non-significant.
The IRRS was calculated for each patient according to the linear combination of these variables weighted by their coefficients of univariate Cox analyses (Table S5). We further split the cases into a high-IRRS subset and a low-IRRS subset. Low-IRRS patients showed significantly better DFS compared with those of high-IRRS (log-rank p < .001) ( Figures 5D-F Figures 5G-I). Additionally, female and more lymph nodes resection implied longer DFS, whereas the advanced N stage indicated the worse in the training cohort ( Figure 5G).
To recognise the immune features responsible for LC recurrence, a heatmap of the 35 candidate variables' immune landscapes and the scatterplot of DFS with corresponding IRRS were plotted ( Figure 6A). We found that helper and cytotoxic T cells, such as CD4+ T cells, CD38+ T cells, CD8+ T cells and CD8+CD133+ T cells, were markedly enriched in the low-IRRS subset, while high-IRRS group showed a higher infiltrating degree of macrophages, such as CD68+ macrophages and CD68+CD163-PDL1-macrophages (both p < .05) ( Figure 6J). Significantly more CD8-CD8 connections were detected in the low-IRRS subset, whereas the high-IRRS subset showed more CD133-CD133 connections and longer distances of CD4-CD20, CD4-CD66b and CD20-CD66b connections (both p < .01) ( Figure 6J). Moreover, high-IRRS patients were in more advanced T stage and N stage compared with low-IRRS populations (both p < .05) ( Figures 6K-Q). Analyses in the testing cohort and entire cohort further validated these findings ( Figure S9-10).
With sensitivity and specificity of 0.655 and 0.724, the IRRS model possessed relatively high predictive accuracy of DFS (AUC 0.738, 95%CI 0.682-0.795) ( Figure 6B). The AUCs for predicting relapse at 1, 3 and 5 years were 0.692, 0.732 and 0.733, respectively ( Figure 6C). The predictive performance of the IRRS model modestly improved when restricted to a specific cTNM stage (Figures 6D-I). The testing cohort (AUC 0.740, 0.707, 0.751 and 0.698 for overall, 1-, 3-and 5-year, respectively) and entire cohort (AUC 0.738, 0.696, 0.738 and 0.722 for overall, 1-, 3-and 5-year, respectively) showed comparable results (Figures S9 and S10). Despite without statistical significance, the IRRS model manifested better predictive ability than the cTNM system S11). Taken together, the IRRS signature presented stable and great performance in separating the recurrence and recurrence-free LC patients after radical resection.

Single-cell transcriptomic landscape of primary tumours
With quality filtering of scRNA-seq data on 4 primary tumour samples, around 33.8 million unique transcripts were obtained from 16,459 cells ( Figure S12). After correcting for reading depth, cells were merged into a dataset with adjustments for batch effects. Dimensionality reduction by UMAP and principal components analysis methods were subsequently conducted. Eleven major cell lineages, including T cells, epithelial cells, natural killer (NK) cells, dendritic cells, macrophages, neutrophils, fibroblasts, B cells, mast cells, stem cells and endothelial cells, were detected based on the expression of classical gene markers ( Figures 7A and B and Table S7). T cells were the The disparities in the clinical characteristics between high and low IRRS subgroups as evaluated by the Chi-square test (L-R). *p < .05; **p < .01; ***p < .001; ****p < .0001; ns, non-significant.

F I G U R E 7
Single-cell RNA sequencing profiling reveals cell composition and cell interaction networks in the lung cancer microenvironment. Eleven major cell lineages as represented in the UMAP plot (A). Canonical gene markers to label clusters in the UMAP most abundant cell types in the primary tumours, followed by epithelial cells and NK cells. Myeloid cells, including dendritic cells, macrophages and neutrophils, were also enriched ( Figure 7C). We further distinguished between tumour cells and non-malignant epithelial cells using several marker genes, including MMP7, MMP9, MMP13 and HOXB2 ( Figure S13). A total of 805 tumour cells and 563 non-malignant epithelial cells were identified eventually.

Spatially adjacent N1-like neutrophils could boost the proliferation and activation of T and B lymphocytes
To further understand the potential biological processes underlying the prognostic effects of neutrophils on spatially adjacent T and B lymphocytes, we first investigated the transcriptomic features of the neutrophil population and three subsets were identified ( Figures 7D and E). The classical subtype was characterised by major expression of canonical neutrophil markers like S100A8 and S100A9. The N1-like sub-cluster displayed features of high expression levels of chemokines and cytokines (e.g., CCL3, CCL4 and C15orf48) and antigen presentation (e.g., HLA-DMB and CD74). The N2-like sub-population was characterised by immunosuppressive (e.g., TGFBI and LAGLS3) and pro-angiogenic (e.g., MMP9) features. [42][43][44] Then the CellChat algorithm was used to investigate the intercellular interaction networks between different clusters of neutrophils and other cell types. First, the overall information flow of each signal pathway was evaluated. The leading intercellular ligandreceptor pairs with neutrophils as signal receivers were distinct among classical (ANNEXIN, ITGB2, ADGRE5) ( Figure 7F), N1-like (MIF, ITGB2, ANNEXIN) ( Figure 7G) and N2-like (APP, MIF, SPP1) clusters ( Figure 7H). N1like cluster showed increased stimulating interactions like MIF-(CD74+CD44), ITGB2-ICAM1 and ADGRE5-CD55, which can promote proliferation and activation of T and B lymphocytes ( Figure 7I). Stimulating interactions with chemotaxis and development functions on dendritic cells and NK cells, such as MIF-(CD74+CD44)/(CD74+CXCR4), also increased in the N1like cluster than N2-like cluster. However, inhibiting interactions like CD22-PTPRC and LAGLS9-CD44/CD45 on T and B lymphocytes increased in N2-like than N1-like cluster ( Figure 7J). Classical neutrophils showed similar interaction networks to N1-like neutrophils ( Figure 7K).
In the external scRNA-seq cohort of LUAD, eleven kinds of cell lineages were identified ( Figure S14A and B and Table S8). Three subsets, including classical, N1-like, and N2-like, of neutrophils were classified based on different expression patterns as well (Figures S14C-F). Major information flow involved in antigen presentation and T cell activation (e.g., MHC-II and MIF) were significantly higher in the classical and N1-like neutrophils ( Figures  S14G and H). In contrast, the N2-like neutrophils showed significantly higher information flow involved in tumour development (e.g., SPP1 and RESISTIN) ( Figure S14I). 45,46 Increased stimulating interactions on T and B lymphocytes, including MIF-(CD74+CD44)/(CD74+CXCR4) and HLA-DRA-CD4, were observed in the N1-like and classical clusters ( Figures S14J and K). However, the N2-like cluster showed increased inhibiting signals like LAGLS9-CD44 ( Figure S14L). Similar results were found in the external scRNA-seq cohort of LUSC ( Figure S15 and Table  S9), validating the above findings.
Then the GSVA analysis comparing N1-like and N2like neutrophils was conducted to further demonstrate the signalling pathways change. A reduction of TGF-β signalling and an increment of interferon (IFN)-α response was observed in N1-like than N2-like cluster, validating that N1-like neutrophils may contribute to stimulating the anti-tumour immune response ( Figure 7L). GO enrichment analysis further illustrated that the modulated genes significantly enriched in inflammatory response and response to IFN-α and IFN-γ biological processes ( Figure 7M). Collectively, an increment of stimulating signals and a reduction of inhibiting signals between N1like neutrophils and major immune effector cells were observed, and the IFN-α and IFN-γ signalling pathways showed a paramount role which drove the functional differences between N1-like and N2-like neutrophils, interpreting the potential biological processes underlying the plot, including CD3D for T cells, NKG7 for natural killer cells, CD163 for macrophages, G0S2 for neutrophils, KRTS for epithelial cells, prognostic effect of neutrophils on spatially adjacent T and B lymphocytes.

Spatially neighbouring M2-like macrophages could blunt the activation of T lymphocytes
The transcriptomic characteristics of the macrophages were profiled, and three subtypes were identified ( Figure  S16A). The pan-macrophage sub-cluster was characterised by markers of cell growth regulation (e.g., EGR1) and tissue-resident (e.g., NR4A3). The M1-like subtype showed major expression levels of IFN (e.g., ISG15) and chemokine (e.g., CXCL10) inducible genes. The M2-like subtype showed similar expression profiles to TAM (e.g., SPP1 and AREG) ( Figure S16B). 47 The overall information with macrophages as signal transmitters were similar among the three subtypes (MIF, MHC-II and APP) (Figures S16C-E). M1-like and M2-like sub-clusters showed an increased signal of angiogenesis than pan-macrophage. Three subsets of macrophages both showed inhibiting signals like LGALS9-CD44/CD45, which may blunt the activation of T and B lymphocytes ( Figures S16F and H). Stimulating signals involved in tumour angiogenesis, like VEGFA-VEGFR2 and VEGFB-VEGFR1, were observed. The M2-like sub-cluster additionally showed increased interactions with tumour cells (e.g., SPP1-CD44) ( Figure S16H). 48 Similar interaction networks were found in the external scRNA-seq cohort of LUAD ( Figure S17) and LUSC ( Figure S18).
We further sought to investigate the signalling pathways change between M2-like and M1-like macrophages in the GSVA analysis, and an increment of genes regulated by MYC and genes involved in angiogenesis was observed ( Figure 7N). GO enrichment analysis further indicated that the altered genes were enriched in chemotaxis and migration of leukocytes and negative regulation of immune effector process ( Figure 7O). Collectively, the above analyses illustrated that spatially neighbouring M2-like macrophages might blunt the activation of T lymphocytes and stimulate tumour growth by promoting angiogenesis.
We herein proposed the model of crosstalk between immune and tumour cells in TME based on the integrated findings from the spatial, receptor-ligand and downstream function analyses (Figure 8).

DISCUSSION
Utilising mIF coupled with deep learning-based spatial assays, we identified various cell phenotypes and mea-sured spatial cell infiltration patterns in the TME of 553 NSCLC cases. We further profiled the transcriptomic features of primary tumours and explored the putative cell interaction networks suggested by spatial analyses. Our findings leveraged insights into the spatial landscape and heterogeneity in TME of NSCLC, underscoring that the spatial relationships between immune and tumour cells might impact clinical outcomes. The correlations between the spatial cell features and prognosis were comprehensively investigated. Proximity among CD20+ B cells, CD4+ T cells and CD38+ T cells was observed, and such spatial configuration has been shown to assist the development of TLS. Moreover, the presence of TLS in NSCLC has been reported to associate with better outcomes. 49,50 Proximity between CD66b+ neutrophils themselves, CD4+ T cells and CD66b+ neutrophils, CD38+ T cells and CD66b+ neutrophils and CD20+ B cells and CD66b+ neutrophils were observed in patients without recurrence than with recurrence, implying the vital role of neutrophils on spatially adjacent cells. To further understand the potential biological processes underlying the prognostic effects of these spatial patterns, then we explored the communication networks mediated by ligand-receptor pairs between different clusters of neutrophils and other cell types and conducted functional analyses using scRNA-seq data. In brief, an increment of stimulating signals and a reduction of inhibiting signals between N1-like neutrophils and major immune effector cells were observed, and the IFN-α and IFN-γ signalling pathways showed a paramount role which drove the functional differences between N1-like and N2-like neutrophils. Consequently, spatial proximity between T and B lymphocytes and N1-like neutrophils, rather than N2-like neutrophils, facilitated the effective anti-tumour immune response, and accordingly, such spatial pattern was associated with a better prognosis ( Figure 8A). Recent studies also showed that neutrophils could stimulate the adaptive immune response by presenting antigens to T cells and producing IFN-γ in early-stage LC. 51,52 They could also coordinate with B cells to mediate the antibody-dependent cellular cytotoxicity process. 53 On the contrary, we found that proximity between Tregs and neutrophils was associated with decreased infiltration levels of CD38+ T cells and CD20+ B cells. We also observed that increased distances between CD66b+ neutrophils and T and B lymphocytes were associated with significantly worse DFS ( Figure 8B). Possible mechanisms may be that the neutrophils (mainly pro-tumour subtype) induced by tumour cells could recruit Tregs to form immunosuppressive TME and promote angiogenesis, contributing to tumour invasion. 54,55 Consequently, our study showed that the prognostic effects of neutrophils could be influenced by their spatial F I G U R E 8 Proposed model of crosstalk between immune and tumour cells of non-small cell lung cancer. Description of patients with a spatially mixed structure among CD4+ T cells, CD20+ B cells and N1-like neutrophils. Spatially adjacent N1-like neutrophils can boost the proliferation and activation of T and B lymphocytes, and such spatial pattern is associated with better disease-free survival (DFS) (A). Depiction of patients with a spatially compartmentalised structure among immune cells (B). More FOXP3+ regulatory T cells and longer distances between neutrophils and CD4+ T cells and CD20+ B cells are observed, and such spatial pattern is associated with shorter DFS. Description of patients with a spatially compartmentalised structure between tumour and CD8+ T cells and have better DFS (C). Greater proximity and strengthened interactions between CD8+ T cells are also found. Depiction of patients with a spatially mixed structure between CD8+ T cells and tumour cells (D). Few CD8+ T cells scatter within these tumours and locate close to CD163+ (i.e., M2-like) macrophages. Spatially neighbouring M2-like macrophages can blunt the activation of T lymphocytes, and such spatial pattern is associated with shorter DFS. location and neighbour cells or milieu, apart from their amounts, whereas few studies have evaluated this issue before. It is acknowledged that due to the lack of distinctive cell surface markers, different functional phenotypes were utilised to define different subsets of neutrophils. 56 Generally, anti-tumour (N1) and pro-tumour (N2) phenotypes of neutrophils were identified and adopted across different cancer types. 57 The phenotypes of neutrophils also alter during the evolution from marrow resident to circulatory and senescent processes. 58 Therefore, the conflicting results reported by different studies could be attributed to different phenotypes of neutrophils. Recently, with the advances of high-resolution approaches, like scRNA-seq, novel and crucial phenotypes of neutrophils were identified, contributing to better understanding their dual roles in TME. 43,59 We found CD133+ CSCs were the most abundant cell populations in TN and closer proximity between CD133+ CSCs predicted significantly worse DFS, similar to the study by Wang et al. 14 that more connections between tumour cells were a poor prognostic factor of LC. Meanwhile, closer distance between CD133+ CSCs was more significantly correlated with recurrence than with more CD133+ CSCs amounts alone, highlighting that spatial analysis could produce a higher-order of useful informa-tion than simply measuring cell amounts. Furthermore, as the second most abundant cell lineages, M2 macrophages were observed to locate close to CD133+ CSCs and CD8+ T cells. Signalling pathways change between M2-like and M1-like macrophages were evaluated, and an increment of genes regulated by MYC and genes involved in angiogenesis was observed. Functional analysis further indicated that the altered genes were enriched in negative regulation of immune effector process, indicating that spatially neighbouring M2 macrophages may blunt the activation of T lymphocytes and stimulate tumour growth by promoting angiogenesis, and accordingly, such spatial pattern was associated with a poor prognosis ( Figure 8D). Previous studies also supported that interactions between M2 macrophages and CD133+ CSCs can mediate the exhaustion of cytotoxic CD8+ T cells and further cause tumour progression. 38,60 Moreover, we found that greater proximity between CD8+ T cells improved DFS. Such trends have been confirmed in other cancer types as well, like colorectal cancer 61 and BC. 62 The findings mentioned above implied that the spatial compartmentalisation between tumour cells and CD8+ T cells ( Figure 8C), instead of mixing ( Figure 8D), was correlated with better DFS. However, the explicit biological events behind the spatial interaction networks were still unclear, and cellular and animal function experimental studies are needed to validate and support our findings.
We then developed a prognostic model integrating spatial and quantitative features that characterised the TME and impacted complex biological pathways. The IRRS was an independent and vicious indicator of recurrence in IA∼IIIB NSCLC patients. The IRRS system displayed relatively high accuracy in predicting recurrence across cTNM stages, contributing to identifying patients at high risk of early relapse and may profit from additional systemic therapies. 63 Moreover, despite without statistical significance, the IRRS model manifested better predictive ability than the cTNM system. The IRRS model presented in the current study has several advantages compared with other prognostic models. First, studies built on cell population estimation algorithms (e.g., CIBERSORT and ssGSEA) only used a limited set of genes. However, gene expression profiles do not necessarily associate with protein levels, and different transcriptional activity among different cell populations can deflate or inflate the appraisal of inactive or active cells potentially. Second, algorithms like CIBER-SORT encountered statistical collinearity when estimating the immune cells. 64 Third, studies based on traditional transcriptome or microarray data could not evaluate the localisation of identified cells.
The immune checkpoint blockade therapies have reformed the treatment of NSCLC over the years. Strikingly, the immune profiles, including immunotherapeutic markers, in different IRRS groups were remarkably different. Patients with low-IRRS displayed an 'immune-hot' subtype, showing the abundant infiltration of CD4+ T cells and CD8+ T cells, and a higher expression of PD-L1 was also detected. As is well known, T lymphocytes are the core of anti-tumour immunity and offer backup sources for immunotherapy. 65 Consequently, with a great abundance of T cells, low-IRRS patients may benefit from immunotherapy potentially. High-IRRS patients were characterised by rich infiltration of macrophages and subsets, and may profit from emerging therapies targeting TAMs. 66 More strikingly, the spatial landscapes between the high-IRRS and the low-IRRS groups were distinct. We found that the crosstalk between CSCs was strengthened, whereas the crosstalk between CD8+ T cells was weakened in high-IRRS patients, hinting at the impaired immune response. Significantly longer distances among CD20+ B cells, CD4+ T cells, and neutrophils were also observed in the high-IRRS group.
Apart from its prognostic value of patients' outcomes, spatial features mining from images also showed potency in early diagnosis of cancer, predicting histologic grade and treatment response. Li and colleagues 67 have established a diagnosis tool for gastric cancer using spatial features from the fluorescence hyperspectral images and reported an overall accuracy of 96.5%. Lagree et al. 68 demonstrated that the spatial characteristics of tumour nuclei in HE images helped improve the classification performance of histologic grade of BC. Xie et al. 69 reported a computational marker derived from HE images that could predict immunotherapy response of LC with an AUC of 0.69. However, most of the models were developed at a single centre and whether these findings could generalise to external datasets remained challenging. Notwithstanding the limitations, previous studies have done a great job of formulating a framework for future studies.
Our study has several distinct advantages. First, the mIF method facilitated detection of multiple protein targets in a single-slide simultaneously, providing a unique perspective of the immune landscape. 3 Additionally, the dataset used in the deep learning model included more than 2700 mIF images of 553 LC patients, that is, a median of 4-5 images were analysed per patient, therefore reducing the bias from blank regions or tissue artefacts that may confound the downstream analyses. Besides, the deep learning model exploited the Delaunay triangulation algorithm without requiring arbitrary definitions of connections, providing a better characterisation of cellular proximity in the TME. 70 Furthermore, the cell communication networks were supported by spatial analyses from mIF graphs and functional analyses from internal and external scRNA-seq datasets.
Despite these advantages, our study has some limitations. First, our study did not include markers for other non-immune components in TME, like fibroblasts. 71 Second, the marker for labelling tumour cells lacked specificity. Although CD133 is well-recognised as a stem cell marker of CSCs, it can also express on normal tissues like immature progenitor cells. 3,72 Also, CD38 is not only expressed on activated T cells, but it can also express on plasma cells and haematological tissues. 73 Likewise, despite being widely adopted in other studies, 74,75 we acknowledged that our method to annotate TN and TS was less objective or precise than directly staining unique markers in TN and TS (e.g., pan keratin and collagen). Third, we did not detect the expression patterns of PD1, which may help estimate the competence of immune cells, like non-exhausted (CD8+PD-1−) or exhausted (CD8+PD-1+) T cells. Moreover, we merely sketched and differentiated the regions of TN and TS in the whole-slide images, without separately analysing the cell composition in the fibrotic or the necrotic regions. Additional limitations included a single institute setting and a lack of external validation of the IRRS model. Future studies in multi-centre settings across populations will be needed to evaluate the strength and validity of our observations. Future studies should focus on the dynamic change in the composition and spatial distribution of cells during the treatment process, for example, neoadjuvant immunotherapy and adjuvant immunotherapy. Furthermore, spatial transcriptomics, in vivo and in vitro experiments are needed for precisely elucidating the rebuilding of the heterogeneous and complex TME in response to different treatment options. As for the IRRS model, future research with treatment information and large sample size is required to assess its value for predicting treatment response.

CONCLUSIONS
We developed a framework to analyse the cell interaction networks in TME based on mIF images, which is efficient for understanding the biological events in which spatial interplays among cells determine functionality. Several spatial architectures, including mixing distributed patterns among CD20+ B cells, CD4+ T cells and N1-like neutrophils, and spatial compartmentalisation between CD133+ CSCs and CD8+ T cells, were underscored to play a vital role in immune response and helped stratify patients who may have an early relapse.

C O N F L I C T O F I N T E R E S T
The authors declare no potential conflicts of interest.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data used to support the findings of this study are available from the corresponding author upon reasonable request.